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^SJ ' An effective wiiite-noise Langevin equation is derived tiiat describes long-time phase dynamics 

of a limit-cycle oscillator subjected to weak stationary colored noise. Effective drift and diffusion 
coefficients are given in terms of the phase sensitivity of the oscillator and the correlation function 
of the noise, and are explicitly calculated for oscillators with sinusoidal phase sensitivity functions 
driven by two typical colored Gaussian processes. The results are verified by numerical simulations 
using several types of stochastic or chaotic noise. The drift and diffusion coefficients of oscillators 
driven by chaotic noise exhibit anomalous dependence on the oscillator frequency, reflecting the 
peculiar power spectrum of the chaotic noise. 
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(~^ . Limit-cycle oscillators are used to model a variety of rhythmic processes in nature. When a limit cycle 

s*/ ' is subjected to noise, the frequency of oscillations changes and the oscillation phase tends to diffuse. 

r^ These effects are quantified by drift and diffusion coefficients and are important in understanding 

•^ the long-time behavior of noisy oscillators. Here, we derive their analytical expressions in terms of 

^ the phase sensitivity function of the limit cycle and the correlation function of the noise, and verify 

'~~' them by numerical simulations using several types of stochastic or chaotic noise. Our formulation will 

. provide a simple and general way to analyze the long-time dynamics of limit-cycle oscillators driven 

^ . by arbitrary weak and smooth colored noise. 

in ; 

[^ ■ I. INTRODUCTION 

^, ; 

lO ' Nonlinear oscillations are ubiquitously observed in nature and various models of limit-cycle oscillators with stable 

^^ , periodic dynamics have been used to describe them [l|, |2|. Some of the many examples are oscillatory chemical 

reactions, cardiac cells, spiking neurons, circadian rhythms, frog calls, passively walking robots, and pedestrians on 

a bridge [l|-[a|- ^ powerful technique for analyzing weakly perturbed limit-cycle oscillators is phase reduction [i1-|J|j 

^ 1 which approximately describes a limit-cycle oscillator possessing multi-dimensional state variables by a single phase 

KJi ' variable. The resulting one-dimensional phase equation is solely specified by the frequency and the phase sensitivity 

'V^ . derived from the original limit-cycle oscillator, which greatly facilitates analytical treatments [3, Q- The diverse 

d ' dynamics of limit-cycle oscillators driven by external forcing or coupled by mutual interactions have been analyzed 

using phase reduction [9l-[l3|. 

Since all systems in nature are subjected to fluctuations, it is essential to incorporate the effect of noise into the 
dynamics of limit-cycle oscillators. It is well documented that noise can induce nontrivial dynamics in oscillator 
systems @tJ, [l3l - [l9| . Synchronization among noninteracting limit-cycle oscillators induced by common or shared 
noisy forcing is a prominent example and has garnered considerable interest in connection with the reproducibility of 
lasers and electronic circuits [20-23], synchrony of spiking neurons [24| - |27| . and large-scale correlated fluctuations in 
ecosystems [28l431j . Phase reduction methods have been extensively used to analyze this phenomenon for limjt-cycle 
oscillators driven by various types of common noise (weak Gaussian noise [32l436j . Poisson random impulses ^37], and 
others 38]). As discussed in |39l - l4l| . phase reduction methods should be applied to noise-driven limit-cycle oscillators 
with care, in order to properly consider the effect of amplitude fluctuations; this is in contrast to the case with smooth 
forcing, where phase reduction can be performed without ambiguity. In [43], a general attempt is made to derive a 
phase equation, which explicitly takes into account the effect of amplitude dynamics, for a wide class of noise. 

In this paper, we analyze noisy limit cycles from an alternative viewpoint, namely, their long-time stochastic 
phase dynamics. In particular, we will derive an effective white-noise Langevin equation that describes the coarse- 
grained phase dynamics of a limit-cycle oscillator driven weakly by sufficiently smooth stationary colored noise. In 
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deterministic systems of limit cycles, e.g., in the synchronization process of coupled oscillators, long-time phase 
dynamics dominate the entire system behavior. Similarly, the long-time behaviors of limit cycles should play crucial 
roles in stochastic oscillator systems and methods for treating them should be developed. Note that the extraction of 
effective dynamics of slow modes has been a classical topic in the theory of stochastic processes (and not just in the 
context of nonlinear oscillators) and various methods such as the use of projection operators and multiscale expansion 
have been developed |43l449l |. 

In the present case, the amplitude effect of the oscillator does not play a significant role at the lowest order 
approximation because its decay is much more rapid than the phase dynamics [41, 42], and hence the conventional 
phase equation holds. We focus on how to obtain drift and diffusion coefficients by specifying the effective Langevin 
equation that gives the long-time phase dynamics of the limit cycle. To this end, we develop a simple theory based 
on the Kramers- Moyal expansion [50, [5l|, which gives the effective drift and diffusion coefficients in terms of the 
phase sensitivity of the limit cycle and the correlation function of the applied noise. Using several types of phase 
sensitivity functions and noisy signals, we demonstrate how the effective drift and diffusion coefficients depend on the 
characteristics of the driving noise. 

II. THEORY 

In this section, we derive an effective Langevin equation describing the long-time phase dynamics of a limit-cycle 
oscillator driven by sufficiently weak and smooth stationary colored noise. We introduce a timescale at which our 
effective description holds, and calculate the drift and diffusion coefficients from the phase sensitivity function of the 
oscillator and the correlation function of the noise. 

A. Model 

We consider a limit-cycle oscillator driven by noise, 

X(i)=F(X)+eC(i), (1) 

where the vector X(t) is the state of the oscillator at time t, F(X) is the intrinsic dynamics of the oscillator, ^(i) is 
the noise, and e is a small parameter representing noise intensity. We assume that Eq. ([T]) has a stable limit-cycle 
solution Xo(i + T) = Xo(i) with period T when the noise is absent (e = 0). The noise is assumed to be smooth, so 
that ordinary rules of differential calculus apply for the variable X(t). More explicitly, we consider the cases in which 
the noise is given by some time-integrated process of (i) stochastic differential equations with Gaussian white noise 
or (ii) ordinary differential equations with chaotic dynamics. 

When the noise is sufficiently small (e ^ 1), the oscillator state can approximately be described using only its 
phase \m^- We first introduce a phase S [0, 27r) on the unperturbed limit-cycle orbit Xo(i) that increases with a 
constant rate (frequency) uj — 2tt/T as 4>(t) = w. This phase (j) can then be extended as a phase field (/'(X) around 
Xo(i) in such a way that (j){t) ~ Vx'?!>(X) • F(X) = uj holds constantly. The dynamics of cj) at the lowest order in e 
are given by 

0(<)=^ + eZ(0(i))e(i), (2) 

where, for simplicity, it is assumed that the noise ^{t) is given only to a single vector component Xi [i = 1, ■ ■ ■ ,N) 
of X, and we denote its intensity by a scalar function £,{t). The 2tt periodic function 
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dX,, 



(3) 

X=Xo(0) 



is called the phase sensitivity [i|-y]j representing a linear response coefficient of the phase <j) to tiny perturbations 
applied to the vector component Xi of the oscillator. Extension to general vector noise is straightforward. 

We assume that ^(i) is a zero-mean stationary random process generated by some noise source, which is smooth, 
temporally correlated, and generally non-Gaussian, with a two-point correlation function C(i), namely, 

(c(i)) = o, (mm)-cit), (4) 

where (• • • ) represents the ensemble average. We further assume that the correlation function decays with a charac- 
teristic time Tc as \C{t)\ = ©(e"!*'/^"^). 



B. Separation of timescales 

Our goal is to derive an effective Langevin equation witli Gaussian wliite noise tliat approximates tlie long-time 
dynamics of Eq. ^. To proceed, we introduce a new slow phase variable by ^(t) — (j){t) — ujt and rewrite Eq. ^ as 

^it) = eZiut + ^jitmt). (5) 

Let Tp — e^^ represent a timescale of the slow phase dynamics of tp, where Tp is much larger than the characteristic 
decay time Tc of the noise correlation C(t), i.e., T^^Tp. 

For sufficiently small e, we can introduce an intermediate timescale t, which is sufficiently longer than the noise 
correlation time, Tc ^ r, but still the slow phase ij) does not change significantly within [t,t + r], namely, 

|V(i + r)-V'(t)|«l. (6) 

This condition implies t <Si Tp = e^^ because |-0(i + t) — ilj{t)\ = 0{€t) for bounded ^(0) and i{t). 
Thus, we have three distinct timescales in our problem, which satisfy 

T^^T <^ Tp. (7) 

The separation of timescales allows us to derive an effective Gaussian white stochastic process from Eq. ([S|) at the 
long timescale Tp that describes the slow dynamics of if) by renormalizing fast fluctuations of the noise ^(i) at the 
short timescale Tc into effective drift and diffusion coefficients. 

C. Effective Langevin equation 

We use a simple Fokker-Planck approximation to the Kramers-Moyal equation ^T] describing the dynamics of a 
probability density function (PDF) P(V', t) of ijj corresponding to Eq. (O, using the periodicity of the phase sensitivity 
function Z{(j)). To this end, we calculate the first- and second-order moments of the slow phase dynamics of -0 during 

[t.t + T], 

Mi(V^,t;T) = (^(i + r)-V(i)), Af2(V,i;r) = ((V^(i + r)-V(i))'), (8) 

where the ensemble average (■ • • ) is taken over noise realizations with fixed "0 and t. Effective drift and diffusion 
coefficients v{il),t) and D(Tp,t), respectively, of the approximate Fokker-Planck equation are obtained from these 
moments at the long timescale Tp{':^ r) by ignoring fast fluctuations of the noise at the short timescale Tc{-^ t). That 
is, we regard r as a small parameter, retain only the 0{t) term, and then formally take the r — > limit; this yields 
the effective drift and diffusion coefficients 

r— »0 T r— ^0 T 

where v{'ip,t) and D{ip,t) will turn out to be constants independent of '0 and t. The resulting approximate Fokker- 
Planck equation 



r^W.') = -— KW'.')1 + -m-2P(*'t) (ID) 



corresponds to an effective Langevin equation, 

V-W ^v + VDf]{t), (11) 

where 'ri{t) is Gaussian white noise satisfying {r]{t)) — and (77(^)77(5)) — S{t — s). Considering the original phase 
variable (f>(t), the effective Langevin equation will be given by 



(t){t)^uj + v + VD7]{t). (12) 

Thus, the effective drift coefficient v gives a noise-induced correction to the raw oscillator frequency uj. The diffusion 
coefficient D gives the effective intensity of the noise. Note that we make no assumption on the oscillator frequency 
u! but assume only that Tc ^ t '^ Tp — e^^ , which can always be satisfied for sufficiently small e. In other words, we 
can always find a scaling region where the above effective Langevin equation is valid, as long as the noise is sufficiently 
weak. 



D. Drift and diffusion coefficients 



To calculate the moments Mi{'ip,t]T) and M2{^,t;T) explicitly, we expand the phase sensitivity function Z[(j}) 
estimated at 4>{ti) = coti + ■ip{ti) as 



Z{ujtl + V'(tl)) = Z{ujtl + ^(t) + 7/'(il) - l/)(t)) 

= Z{ujt, + ^{t)) + Z'iiut, + ijit))mh) - ^{t)} + Oimh) - 7/;(i)}2), 
where Z'{(j)) — dZ{<f>)/d(f), and we integrate Eq. ([5]) as 



(13) 



t + T 



^{t + r) - ij{t) = e / dtiZ{ujti + iP{t))^{ti) 



t + T 



dhZ'{ujh + i>{t))mti) ~ ^WK(ii) + 0{e{i,{h) - i:{t)}). 



(14) 



By iterative substitution, namely, by inserting ipiti) — ip^t) — e J^^ dt2Z'{ujt2 + ip{t))^{t2) + 0{e^) obtained from the 
above equation into its second term, we obtain 

l-t + T 

ipit + r) - ij{t) = e / dtiZ{ujti + '0(i))C(ii) 

rt+T rti 

dh / dt2Z'iujh + ^{t))Z{u:t2 + ^{m{ti)^{t2) + 0(e', T^), 



(15) 

where we have used the fact that \'(p{t + t) — ipi^)] = 0(er). Taking the ensemble average of this expression over the 
noise, the moments Mx{ij),t\T) and M2(ip,t;T) can be calculated up to 0{e^,T^) as 



Mi(V',i;T) 



M2(V',i;r)=e2 



t + T pti 

dti dt2Z'{ujti+ij)Z{ujt2+ij)C{t~ti) 
t Jt 



t + T i^t + T 

dh dt2Z{ujti+ip)Z{ujt2 + tp)C{t-ti) 






Now we take r as an integer multiple of T = 27t/uj, namely, 

T = nT ~ 2n7r/w 



(16) 
(17) 

(18) 



with some appropriate integer n (= 1, 2, • • • ) such that Tc -€1 t = nT <ti Tp = e "'^ is satisfied (we may simply set n = 1 
if Tc <^ T). Using the periodicity of the phase sensitivity function, the moments can be written as 



Afi(V',i;T)=eV 



M2(^,i;r)=eV 



— / dsC{s) / d9Z'{e)Z{e-ujs) 
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dsC{s) / dez{e)z{e-ujs) 
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+ 0(e^r2), 



(19) 
(20) 



both of which turn out to be constants (see Appendix for calculations) . 
From Eq. (J9]), the effective drift and diffusion coefficients are obtained as 



— / dsC{s) / dez' 



{e)z{e-ujs) 



and 



D = e' 



1 

2^ 



dsC{s) / deZ{9)Z{9-us) 
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(21) 



(22) 



The colored noise gives constant contributions of O(e^) to both v and D. In [53j and |54| . similar weak-noise expansion 
methods for the phase dynamics of noise-driven limit cycle oscillators are used to estimate their Lyapunov exponent 
(Eq. (8) in [53J, which generalizes the result in [33| for Gaussian noise) or the variance of periods. 



E. Fourier representation 

The effective drift and diffusion coefficients can be expressed concisely using the Fourier representation of the phase 
sensitivity function, 



z{e)^ Y, Zi^''\ (23) 



as well as the power spectrum of the noise ^(i), 

/oo 
City^Ht. (24) 

-00 

Because ^(t) is stationary, the correlation function satisfies C(t) — C{~t)^ so that the power spectrum I{^) can be 
expressed as 

I{n) = 2 / C{t) cos{nt)dt = 2Rex(fi), (25) 

Jq 

where 



00 



Xm - / C(i)e*"*di (26) 

Jo 

is the Fourier-Laplace transform of the correlation function. Since x(r2) is analytic in the upper-half of the complex 
plane (Re i7 > 0), we can express its imaginary part using the Kramers-Kronig relation as 

Imxm = -P.V. r ^^dz = \m, (27) 

where 

1 f^ I(z) 

I{n) ^ -RV. / TT^^dz (28) 

n J_^il-z 

is a Hilbert transform of the power spectrum I{Vt) (see e.g., [53|)- Inserting these equations into Eqs. (PT|) and (|22p . 
the drift and diffusion coefficients v and D can be expressed as 

2 '^^ 
^=Y E (*^)l^£p{n^^)+^/(^^)}, (29) 

t=-oo 

and 

CXD 

D^t^ Y. \Zi\^H^^)- (30) 

f=-oo 

Thus, V and D can be calculated from the power spectrum I{fl) and its Hilbert transform I{^). This is convenient 
because the phase sensitivity Z((j)) often contains only lower harmonic components. 

III. EXAMPLES 

In this section, we numerically verify the accuracy of the effective white-noise phase Langevin equation (jlip for 
several types of colored noise generated by stochastic processes and deterministic chaotic systems. We compare 
the effective drift and diffusion coefficients given in Eqs. (|2ip and ([2^ . respectively, with those obtained by direct 
numerical simulations of the original phase model, Eq. ([5]). 
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FIG. 1: Type-I sinusoidal phase sensitivity Zi{(j)), Type-II sinusoidal phase sensitivity Zii{<j)), and Type-E phase sensitivity 
ZEifj)) obtained from a Morris-Lecar spiking neuron model. 



A. Phase sensitivity functions 

We consider the following examples of phase sensitivity functions (see Fig. [1]) : 
1. Type-I sinusoidal function with only a positive lobe, corresponding to limit cycles near saddle-node bifurcation [3, 



1 ype-i s 

mm, 



zm 



(31) 



2. Type-II sinusoidal function with positive and negative lobes, corresponding to limit cycles near Hopf bifurca- 
tion ^], 



Zii[<i>) 



(32) 



3. Phase sensitivity function Ze{<J)) of the Morris-Lecar neuron model near homoclinic bifurcation (see Appendix), 
where the function is not sirnply sinusoidal but contains higher-order harmonics. It can be calculated numerically 
by the adjoint method [111, l55|. We refer to this function as Type-E. 

The first two functions Ziji{(j)) are generic in the sense that they can be derived analytically from the normal forms of 
limit-cycle oscillators near the respective bifurcation points by appropriate coordinate transformations |12l . |55| . The 
third function Ze{4>) is model-dependent, but is a typical example of the phase sensitivity near a homoclinic bifurcation 
point. ZE{(f>) tends to be dominated by an exponentially decaying part resulting from the linear dynamics near a 
saddle point as the bifurcation point is approached, and therefore a simple exponential function with a discontinuity 
is proposed as a generic form of the phase sensitivity in [55|. We do not however use this form to avoid unnatural 
effects of the artificial discontinuity. 



B. Measuring the coefficients 

We estimate the effective drift and diffusion coefficients v and D by direct numerical simulations of Eq. ([5|) and 
compare them with the respective theoretical values, Eqs. (PT|) and (|22p . The solution to the effective Fokker-Planck 
equation pH)) from a delta-peaked initial condition P{ip,0) — 5{^Jj) is simply a Gaussian wave packet. 
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FIG. 2: (a): Sample paths of the slow phase i/)(f) driven by the Ornstein-Uhlenbeck noise (100 realizations). Oscillator 
frequency lj — 1 and the phase sensitivity function is Zii{<j>). Noise correlation time Tc — 1 and noise intensity e = 0.1. Solid 
lines represent the mean value and broken curves indicate the mean value ± the standard deviation, (b) and (c): Temporal 
growth of the mean (b) and the variance (c) of the phase 7/)(i) averaged over 200, 000 realizations for a; = 1, 2, 3, and 4. The 
other parameters are the same as in (a). 



whose moments are given by 



m))^vt, {mt)~m))r 



Dt. 



(34) 



Thus, we can measure v and D from slopes of the mean and the variance of the phase plotted as functions of t. 

For example. Fig. [2ja) displays typical sample paths of Eq. ([5]) with the Type-II function Z// (</>). The evolution 
of the slow phase ?A(i) = (t){t) — tut is plotted for 100 realizations of the OU noise (explained below). The broken 
line represents the mean path averaged over 200, 000 realizations, which shows negative drift induced by the finite 
correlation time Tc = 1 of the noise. Figures Hfb) and (c) display the mean and the variance of the oscillator phase, 
respectively, averaged over 200,000 realizations for differing values of a; and for the Type-II Z((/)), all of which clearly 
show linear dependence on time t, whose slopes yield v and D. 



Ornstein-Uhlenbeck noise 



We first consider the case in which the colored noise ^(i) obeys the OU process, 

Tc Tc 



(35) 



where rjit) is zero-mean Gaussian white noise whose correlation function is given by {'i]{t)ri{s)) = 5{t — s). This OU 
process generates colored Gaussian noise ^(i) with a stationary PDF 



1/2 



PiO - (^) eM-Tce) 

and an exponentially decaying correlation function 



2Tr \ Tc 



(36) 



(37) 



Thus, the characteristic decay time of the noise correlation is Tc- In the limit Tc — > 0, C(i) converges to the Dirac 
delta function S{t), so that ^(i) converges to Gaussian white noise of unit intensity. The power spectrum of ^{t) and 
its Hilbert transform are given by 
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(38) 
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FIG. 3: Effective drift and diffusion coefficients v and D, respectively, piotted against osciilator frequency uj [(a) and (b)] and 
correlation time Tc [(c) and (d)] for Type-I, II, and E phase sensitivity functions Zi{(j)), Zn{(j>), and Zsifj}), respectively, driven 
by Ornstein-Uhlenbeck noise. Noise intensity e = 0.1 in all cases. Noise correlation time Tc = 1.0 in (a) and (b), and oscillator 
frequency to = 1.0 in (c) and (d). Data obtained by direct numerical simulations are compared with theoretical curves. 



By inserting Eq. (|38)) into Eqs. (|2T|) and ((22)) . the drift and diffusion coefficients are expressed as 
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+ (WT,^)2 ' 



and 
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Note that v is always non-positive and vanishes in the white-noise limit [tc — ^ 0). That is, the OU noise <^(i) always 
tends to slow down the oscillator for arbitrary (smooth) phase sensitivity functions even if (f (t)) = holds on average, 
which agrees with the result previously obtained by Galan [52] (however, [521 ^^es a different definition of the OU 
process). Similarly, the diffusion coefficient D is maximized in the white- noise limit. 
For the Type-I phase sensitivity function Zj{(f)), v and D are explicitly calculated as 



4 1 + (LJTcy 
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(41) 



and for the Type-II Zjj{(j)) as 



D^^ , , , ,, . (42) 



4 1 + (cjr,)2 2 1 + (cjTc) 



Note that v is the same for both Zi{(j)) and Zii{(t)), whereas D for Zi{(j)) is larger than that for Zn {(/)). This can 
easily be seen from the Fourier representations; the only difference between Zi{(j)) and Zij{(j)) is that Zi{(p) has a 
non-vanishing constant component Zq = 1/2. For the Type-E phase sensitivity Ze{4'), we numerically integrate 
Eqs. dill) and ^^ to obtain v and £>. 

In numerical simulations, the noise correlation time is fixed at Tj, = 1 and the noise intensity at e = 0.1. Figures[3l^a) 
and (b) plot v and D as functions of the oscillator frequency uj for the three types of phase sensitivity (averaged over 
200,000 realizations) and compare them with theoretical values, indicating good agreement. The drift coefficient v for 
Zi{(j)) and Zii{(j)) coincide with each other and are minimized at oj = t~^ = 1. The diffusion coefficient D for Zj{(j)) 
and Zii{(f>) differ from each other and decrease monotonically with Tc- In particular, D for Zii{(j)) (more generally 
for Z{(l)) without a constant component Zq) tends to vanish at large w, indicating that the long-time phase diffusion 
of Typc-II oscillators can be very small when the oscillator frequency is large. The numerical values and theoretical 
values of v and D for Ze{(I)) are also in good agreement. 

D. Noise generated by a damped noisy harmonic oscillator 

Next, we consider colored noise ^(t) generated by a damped noisy harmonic oscillator (hereafter referred to as 
DNHO noise), 

x = uJoy-'yx + 'yr]x{t), y ^ -ujqx - jy + jr]y{t), (43) 

where ?7x(0 and riy{t) are mutually independent Gaussian white noise satisfying {rix{t)rix{s)) = {Vy{'t)''ly{s)) = S{t — s) 
and {r]x{t)riy{s)) — 0. The parameter loq is the frequency of the harmonic oscillations and 7 is the damping constant. 
This process yields two-component noise with a Gaussian stationary PDF, 



P{x, y) = — exp 

7r7 

and a correlation function C{t) of x{t) with oscillatory decay. 
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(44) 



C{t) = {x{t)x[Q)) = ^e-'^l*l cos(a;oi). (45) 

Thus, the correlation time is given by Tc — ^^^ ■ We use this x{t) as the noise ^(i) given to the oscillator. The power 
spectrum of x{t) and its Hilbert transform are respectively given by 

^(^) = V f ^-777— ^ + ^-777 ^^ ' (46) 

2 V7^ + (S2 + wo) 7^ + (S2-a;o) / 

and 

Effective drift and diffusion coefficients v and D can be analytically calculated for the Type-I phase sensitivity 
Zi{(f)) as 

6^7 / LU — LOq UJ + LUf) 



7^ + (w — wo)2 72 + (cj + wo)^ 
^ " ^ U2 + (t^-t^o)2 ^ ^^T^ ^ 72 + (c. + L.o)2 ' ' ^^^^ 
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FIG. 4: Effective drift and diffusion coefficients, v and D respectiveiy, piotted against osciliator frequency u and correlation 
time Tc for phase oscillators with Type-I, II, and E phase sensitivity functions Zi{(j>), Zii{(j>), and Ze{<I>), respectively, driven 
by the damped noisy harmonic oscillator noise. Noise intensitye = 0.1 and noise frequency ujq = 2.0 in all cases. Tc = 1.0 in (a) 
and (b), and ui = 1.0 in (c) and (d). Data obtained by direct numerical simulations are compared with theoretical curves. 



and for the Type-II phase sensitivity Z2{(j)) as 



L> = 



6^7 


/ 


<jj — 


Wo 


+ 




W + Wq 


8 


\y 


+ (w 


-^o)2 


r 


+ (w+Wo)2 


£2^2 


f 


1 




+ 




1 ^ 


4 


v7^ 


+ (cj- 


-^o)2 


7^ 


+ {uj + wo)2y 



(49) 



In the cjg — > limit, DNHO noise returns to the OU noise, so that v and D converge to the corresponding results 
for the OU noise. Note that values of v coincide again between Zi{(j)) and Zii{(l)), vifhereas those of D differ between 
the two cases. Values of v and D for the Type-E phase sensitivity ZEifj)) ^-re calculated by numerically integrating 
Eqs.dUl) and (|22)) . 

Figure [4] plots v and D obtained by direct numerical simulations of Eq. ([5]) (averaged over 200, 000 realizations), and 
the data are compared with the theoretical results. The parameters e — 0.1 and wq = 2 are fixed and the oscillator 
frequency ui or the noise correlation time Tc — 7~^ is varied. In Figs. Slja) and (b) their dependence on uj with fixed 
Tc = 1 is shown. In contrast to the OU case, v can take positive and negative values for all Z{(f)). D does not 
monotonically decrease but exhibits a peak (at tu = luq for Zj{(j)) and Zjj{(j))) implying some type of resonance effect. 
D for Zj{(j)) is again larger than that for Zii{(j>). Figures HKc) and (d) show the dependence of v and D on the noise 
correlation time Tc with fixed oscillator frequency w = 1. The drift coefficient v can take positive values for Zj((j)) 
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and Zij{(f>). D decreases monotonically for all types of the phase sensitivity. In all cases, numerical and theoretical 
results are in agreement. 

E. Noise generated by a chaotic Lorenz model 



The Lorenz model [56 



p{—x + y), ij = —xz + qx — y, z = xy — rz (50) 



generates a typical chaotic time sequence. We use parameter values p — 10, q ~ 28, and r — 0.9 and apply the 
normalized time sequence of x{t), 

m - ^-^^. (51) 

to the phase model as the colored noise ^(t), where (• • •) denotes the long-time average. Figures [5fa) and (b) show 
the correlation function and power spectrum of the noise, respectively. The correlation function exhibits oscillatory 
decay with several characteristic frequencies, which appear in the power spectrum of the noise as sharp peaks. 

For simplicity, we consider only the Type-II phase sensitivity Zjj{(jj). Figure [5jc) and (d) compare the drift and 
diffusion coefficients v and D obtained by direct numerical simulations of Eq. ([S]) (averaged over 100,000 realizations) 
with the Lorenz model and the theoretical values calculated from the correlation function C(i), which are in agreement. 
The noise intensity is fixed at e = 0.1 and the oscillator frequency w is varied. The drift and diffusion coefficients v 
and D show interesting peculiar dependence on uj. As w increases, v increases rapidly and then suddenly decreases, 
and this is repeated several times. D exhibits a few sharp peaks, indicating that phase diffusion due to the Lorenz 
noise can be strongly enhanced at some particular values of the frequency uj. 

These results, in particular the behavior of D, can easily be understood from the Fourier representation, Eq. (|30l) . 
Since Zjj{(j)) has only the first harmonic component, Z±i = ±i/2, Eq. ([5D|) simply gives D = e^[I{uj) + /(— w)]/4 = 
e^/(w)/2, namely, D is simply proportional to the power spectrum itself. In fact, we can see that the curves in 
Figs.[5jc) and (d) are identical except for the scaling factor e^/2. Moreover, from Eq. ((29)) . we can see that the sudden 
rise and fall of v is due to the Hilbert transform I{fl) of the power spectrum I{fl) near its sharp peaks, which gives a 
contribution l/(fl — fl') if the peak is approximated by a Dirac 6 function d{fl — fl'). 



The Rossler oscillator [56 



F. Noise generated by a chaotic Rossler oscillator 



-y — z, y = —X + ay, z = b + xz — cz (52) 



is another typical example of low-dimensional chaos. We fix the parameter values at a = 0.3, b — 0.2, and c — 5.7, 
where the Rossler oscillator possesses a "funnel" attractor. The cc-component is normalized as in Eq. (jSTj) and applied 
to Eq. ^ as the noise £,{t). Figures UJ a) and (b) show the correlation function and power spectrum, exhibiting 
oscillatory decay and sharp peaks similar to the Lorenz model. 

We consider only the Type-II phase sensitivity Zjj((j)) again. Figure [SJc) and (d) compare the drift and diffusion 
coefficients v and D obtained by direct numerical simulations of Eq. ([5]) (averaged over 100,000 realizations) with 
those calculated from the correlation function. Noise intensity e = 0.05 is fixed and oscillator frequency lu is varied. 
Similarly to the case of the Lorenz model, v and D indicate interesting complex dependence on the oscillator frequency 
w, reflecting the peculiar power spectrum of the Rossler oscillator. There is again a good agreement between the 
numerical and theoretical values. 

IV. SUMMARY 

We derived an effective white-noise Langevin equation that describes the long-time phase dynamics of a limit-cycle 
oscillator driven by general non-Gaussian colored noise. Effective drift and diffusion coefficients were calculated from 
the phase sensitivity of the oscillator and the correlation function of the noise. The results were verified using several 
types of colored noise sources, i.e., the Ornstein-Uhlenbeck process, the damped noisy harmonic oscillator, and the 
chaotic Lorenz and Rossler models. 
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FIG. 5: (a) Correlation function C{t) and (b) power spectrum I{uj) of the normalized x variable of the Lorenz model, (c) Drift 
and (d) diffusion coefficients of the phase with Type-II phase sensitivity Zu (</>) driven by the Lorenz model. 



Our analysis gave general expressions for drift and diffusion coefHcients, applicable to general limit-cycle oscillators 
driven by arbitrary weak smooth noise. In a previous study [52| , Galan calculated the frequency shifts of limit-cycle 
oscillators (the effective drift coefficient v in our notation) driven by colored Ornstein-Uhlenbeck noise, and pointed 
out that the frequency shift is always negative for arbitrary phase sensitivity. In contrast, for other types of noise, 
frequency shifts can also be positive, so that the noise may increase the frequency of the driven oscillator. We can 
also calculate the effective diffusion coefficient D, which directly reflects the power spectrum of the driving noise. In 
particular, for chaotic noises, D exhibited sharp peaks, indicating that the phase diffusion can be greatly enhanced 
for peculiar frequencies of the limit-cycle oscillator. 

The effective white-noise Langevin description enables us to use the powerful classical methods for stochastic pro- 
cesses ,50, 51] and thus provides a general framework for analyzing the long-time behavior of limit-cycle oscillators 
subjected to noise. Important future topics will include generalization of the present results to multi-dimensional situ- 
ations and incorporation of deterministic external forcing (e.g., periodic) and mutual interactions. It is expected that 
the combined effect of colored noise and other external perturbations or mutual interactions may lead to qualitatively 
new dynamics. 
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FIG. 6: (a) Correlation function C{t) and (b) power spectrum /(cj) of the normahzed x variable of the Rossler oscillator, (c) 
Drift and (d) diffusion coefficients of a phase oscillator with Type-II phase sensitivity Z//((/)) driven by the Rossler model. 



V. APPENDIX 

A. Derivation of Eq. ([T9|l from Eq. ([T6|l 

We set r — nT with n being an integer. The right-hand side of Eq. (|16p can be rewritten as 

dh / dt2Z'{ujti + '4}{t))Z{ujt2 + ip{t))C{h - 12) 
t Jt 

t+T nti-t 

dti / dsZ'{ujti+tp{t))Z{uj{ti~s)+tp{t))C{s) 
t Jo 

Pt + T 

dsC{s) / dtiZ'{ujti + i}}{t))Z{uj{ti -s) + ^P{t)) 
Jt 

27r 



27r 



dsC{s) / d9Z'{0)Z{0-ujs), 



(53) 
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where we have used 

ft + T 



dtiZ'iojti + -0(i))Z(a;(ti - s) + ipit)) 
t 

71 — 1 pX 



V / dtiZ'{Lj{ti+t + jT) + i^(t))Z{Lj{ti-s + t + jT) + tP{t)) 

= V — / d0Z'{e + ut + 2ttj + tlj{t))Z{9 - ujs + ujt + 2ttjT) + jp{t)) 

= — / dez'{e)z{e-ujs) = — dez'(e)z(e~ujs). 

27r Jo 27r Jq 



(54) 



Substituting these results into Eq. (IT51) yields Eq. (|19l) . 



B. Derivation of Eq. ([20]) from Eq. (fTf)) 

Setting T = nT, the right-hand side of Eq. (IT71) can be transformed as 

dti / dt2Z{ujti + i:{t))Z(ujt2 + ^-(0)^(^1 - ^2) 
t Jt 

t+T ntt-t 

dti / dsZ{ujti + iJj{t))Z{u(ti -s)+ 'ip{t))C{s) 

00 pt-\-T 

dsC{s) I dtiZ{ujti + i;{t))Z{uj{ti - s) + 2P{t)) 

-00 -/i 

00 p2tt 



dsC(s) / d0Z{0)Z{0-ujs) 
-00 Jo 



(55) 



where we have approximated the range of the integral of the correlation function C{s) over [ti — t — t. ti — t] as 
[—00, +00] by assuming that the decay time of C{s) is much shorter than t and that ti — t — t < and ti — t > 
hold. In deriving the final expression, we used 



t+r 

dtiZ{u}ti + ip{t))Z{uj{ti - s) + tpit)) 
t 

n — 1 f,T 



= V / dtiZ{uj{ti+t + jT) + 'ijj{t))Z{uj{h-s + t + jT) + 'ilj{t)) 

n-1 rp p2Tr 

= V — / d0Z{0 + u}t + 2Trj + i'{t))Z{0 -ujs + ujt + 2ttj + tp{t)) 

^ — ' 2tt If, 

= y— / d0Z{0)Z{0-iJs)^-- d0Z{0)Z{0-ujs) 

^ 27r Jo 27r Jo 

= _ / d0Z{0)Z{0-LJs). (56) 

27r Jo 



/o 
We obtain Eq. (PH)) by substituting the above result into Eq. (IT71) 
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The Morris-Lecar model 



The Morris-Lecar model of a spiking neuron is given by the foUowing set of two-variable ordinary differential 
equations |5ll9l.l55|: 



CV{t)=gcam^(y){Vca 



V) + gK{VK -V)+ gUVL -V) + I, 



w(t) = (j) 



where 



T^(V^) 






(57) 



1 + tanh 
1 + tanh 



T^{V) 



cosh 



2Va 



V -Vi 

Vi 
-I -1 



(58) 



Here, V represents membrane potential and w is an activation variable for potassium. The parameter values are chosen 
as = 0.23, QL = 2.0, gca = 4.0, gK = 8.0, C = 20.0, Vk = -84.0, 14 = -60.0, Vca = 120.0, V^ = -1.2, V^ = 18.0, 
V3 — 12.0, V4 = 17.4, and / = 37.5 f9]. This model exhibits limit-cycle oscillations via homoclinic bifurcation near 
/~35. 

We set the origin of the phase (f) — aX the point where V exceeds from below. The phase sensitivity function Ze{(J)) 
for this model can be numerically obtained by the adjoint method as explained in [ll|, [S^]- It has an exponentially 
decaying part, which tends to dominate the whole function as the parameter I approaches the bifurcation point. The 
above set of parameter values gives a fixed frequency w = 0.198. However, note that we may still set a; arbitrarily as 
in Figs. [3] and m by rescaling the time appropriately {Ze{4>) is not affected by time rescaling). 
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